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Abstract 

A new procedure aiming at folding a powder diffraction 2-D into a 1-D scan 
is presented. The technique consists of three steps: tracking the beam centre 
by means of a Simulated Annealing (SA) of the diffraction rings along the 
same axis, detector tilt and rotation determination by a Hankel Lanczos Sin- 
gular Value Decomposition (HLSVD) and intensity integration by an adaptive 
binning algorithm. The X-ray powder diffraction (XRPD) intensity profile of 
the standard NIST Si 640c sample is used to test the performances. Results 
show the robustness of the method and its capability of efficiently tagging the 
pixels in a 2-D readout system by matching the ideal geometry of the detector 
to the real beam-sample-detector frame. The whole technique turns out in a 
versatile and user-friendly tool for the 2?? scanning of 2-D XRPD profiles. 

* Corresponding author, E_mail: massimo.ladisa@ic.cnr.it, Phone: +39 0805929166, Fax: 

+39 0805929170 



I. INTRODUCTION 



X-ray powder diffraction (XRPD) technique is nowadays a well known tool to study 
crystalline properties, which provide important information for applications in many fields 
such as crystal structure solution, microstructure characterization, phase quantification and, 
recently, also in nanotechnology. In these fields, charge-coupled device (CCD) and imaging 
plate (IP) detectors have been widely used for collecting X-ray diffraction images, expecially 
because of their fast readout system (Hanley et al. 2005; Hammersley et al. 1997, Muchmore 
1999, Phillips et al. 2002, Suzuki et al. 1999). When the analysis of the 1-D powder pattern 
is rather sofisticated, as for the Rietveld techniques, a reliable pre-processing aimed at 
enhancing the quality of XRPD data is required. When working with 2-D detectors, like a 
CCD or an IP system, a crucial step is the folding of the 2-D powder diffraction collected 
image into a 1-D scan. In the present paper a new procedure aiming at folding a powder 
diffraction 2-D into a 1-D scan is presented. 



II. THE METHOD 

The ideal X-ray powder diffraction geometry is cylindrical along the beam-sample axis. 
The diffraction cone with an opening angle 2i? and the apex on the sample reaches the 
imaging plate (a CCD detector for instance). Ideally the imaging plate would be orthogonal 
to the beam-sample axis and it would collect a series of circles with radii corresponding to 
different momentum transfer: the larger the radius the higher the momentum transfer (i.e. 
the resolution). Unfortunately the imaging plate axis rarely fits the ideal geometry and 
the circles become ellipses with an eccentricity e and a rotation 0o parametrizing tilt and 
rotation of the detector frame. Furthermore the beam rarely reaches the geometrical centre 
of the plate and a global shift of the rings is needed to perform a correct 2d scanning. 
The problem reduces to the nonlinear fit of the parameters appearing in the cartesian equa- 
tion of an ellipse: x 2 / a 2 + y 2 /b 2 = 1, where a/b is the major/minor semiaxis. For a generic 
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beam-sample-plate geometry the following equations hold: 



x = x c + Ax c + a cos(0 — 0o ) / \A + e sin 2 (0 — O ) , 
V = Vc + Ay c + a sin(0 - o )/\/l + e sin 2 (0 - O ) , 



(1) 



being Ax c , Ay c the discrepancies along x,y axes between the beam centre and the imaging 
plate one (x c , y c ). e is the ellipse eccentricity a 2 /b 2 — 1. The cartesian equation for the ellipse 
finally reads: 



1 + eS 2 -l = 2p 



Ax c C Ay c S 



b 2 



Ax 2 c Ay 2 c 



b 2 



(2) 



where p = \J (x - x c ) 2 + (y - y c ) 2 , S = sin(0 - </> ), C = cos(0 - <f> ), and = O + 
arctan ((y - y c )/ (x - x c )). 

In Eq. (2) the l.h.s. is 0(e) while the r.h.s. is O ( 



and, since usually 



e \/ ^ — -, Eq. (2) finally decouples. It turns out that the problem of finding the 

V a b 

rings centroid can be solved independently from the e/0 o determination which is indeed a 
fine tuning with respect to the former task. Needless to say that the attempt of finding the 
rings centroid cannot be accomplished with by using the eq. (2) itself where both e and 4>o 
are unknown at this stage. 

The decoupling of the two problems, i.e. centroid hunting and ellipses fitting, is the main 
advantage of our method in comparison to the ones nowadays available in the literature 
(see Hammersley et ai, 1996) where the full set of parameters of an equation analogous to 
Eq. (2) is computed. 



A. Beam centroid hunting 

This task is carried out by using a Simulated Annealing (SA) algorithm (see Metropolis 
et al. 1958, Pincus 1970, Kirkpatrick et al. 1982, Kirkpatrick et al. 1983, Kirkpatrick 
1984) to search for a minimum of the XRPD rings centroid position on the imaging plate. 
As known, the algorithm employs a random search which not only accepts changes that 
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decrease objective function /, but also some changes that increase it. Changes decreasing 
/ are accepted with a probability p = exp(— Sf/T), where 5f is the decrease in / and T is 
a control parameter, which by analogy with the original application is known as the system 
temperature irrespective of the objective function involved. Here the objective function / 
is the sharpness (defined e.g. as the inverse normalized variance) of the radial distribution 
function of the highest intensity rings recorded on the imaging plate and the system tem- 
perature is the domain size around the guess beam centre during the random search. / is 
computed with respect to the guess rings centre. In Fig. 1 the sharpness of the radial distri- 
bution function is shown before (top right) and after (bottom left) the SA. The increasing 
of / sharpness is dramatic and it follows the correct positioning of the beam centre onto the 
imaging plate (see the two spots on the top left image of Fig. 1). 



B. Detector tilt and rotation determination 

Once the beam centre has been positioned, i.e. Ax c = Ay c = 0, the eccentricity e 
and rotation O computation is straightforward; infact the eq. (2) reads 1/p 2 ~ l/a 2 x 
1 + e sin 2 (0 — O ) . The detector tilt and rotation are determined by using a subspace- 
based parameter estimation method called Hankel Lanczos Singular Value Decomposition 
(HLSVD) technique (Laudadio et al, 2002). The HLSVD method works as follows. Let us 
model the 1/p 2 by samples collected at angles </>„, n = 0, ... N — 1 as the sum of K 
exponentially damped complex sinusoids 

K 

1 /Pl-Yl a k exp(-4 <Pn) C0S[27T f k (j) n + tp k \ , (3) 

k=i 

au is the amplitude, ipk the phase, dk the damping factor and fk the frequency of the k th 
sinusoid, k = 1,...,K, with K the number of damped sinusoids. This multiparametric 
functions family has been successfully used (Laudadio et a/., 2002) for modelling NMR 
signals, because these functions are intrinsically related to the physical description of a 
NMR signal. It has been shown, however, that this family has a great approximating 
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power for many kinds of signals, including some of the natural functions used to describe 
crystal diffraction patterns, namely the Dirichlet kernel, related to the the Laue interference 
function (see Beylkin & Monzon 2005, Sec. 5.1). Practical tests (Ladisa et al. 2005a; 
Ladisa et al. 2005b; Cervellino et al. 2005a; Cervellino et al. 2005b) have also confirmed 
their effectiveness. To estimate the parameters involved in the approximation we proceed 
as follows. The N data points defined in (3) are arranged (see Beylkin & Monzon 2005, 
Sec. 3) into a Hankel matrix H d = Hf x l M = l/p 2 m+l , m = 0, . . . , M - 1, I = 0, . . . , L - 1 
, with L + M = N + 1 (M~L~ N/2). The spectral decomposition of the Hankel 
matrix is performed by SVD. SVD gives this method its main advantage, that is flexibility, 
since the number of parameters is not fixed and it is determined only by the desidered 
agreement level between the model and the real peak profile shape. SVD of the Hankel 
matrix is the decomposition H LxM = Ulxl^lxmVmxMi where E = diag{Ai, A 2 , . . . , A r }, 
Ai > A2 > . . . A r , r = min(L, M), U and V are orthogonal matrices and the superscript H 
denotes the Hermitian conjugate. A fast SVD decomposition is achieved by the Lanczos 
bidiagonalization algorithm with partial reorthogonalization. This algorithm, based on FFT, 
computes the two matrix-vector products which are performed at each step of the Lanczos 
procedure in 0((C + .M)log e (£ + M)) rather than in 0{CM). In order to obtain the 
signal subspace, the matrix H is truncated to a matrix Hk of rank K Hk = Uk^kV^ , 
where U K) V K) and Y^ K are defined by taking the first K columns of U and V, and the 
K x K upper-left matrix of S, respectively. As subsequent step, the least-squares solution 
of the following over-determined set of equations is computed E H ~ y^ ottom ) ^ w here 
yjbottom) an( j y(t°p) are derived from V K by deleting its first and last row, respectively. The 
K eigenvalues Zk of matrix E are used to estimate the frequencies fk and damping factors 
dk of the model damped sinusoids from the relationship Zk = exp dk + i2iifkj A$ , with 
k — 1, . . . , K. Values so obtained are inserted into the model equation (3) which yields the 
set of equations 

K 

l /pl-J2 a k exp(-4 <Pn) cos[27r/ fc 0„ + tp k \ , (4) 
fc=i 



with n — 0, . . . , N — 1. The least-squares solution of (4) provides the amplitude dk and 
phase 0k estimates of the model sinusoids which are used in the next step. 
The eq. (4) and the l.h.s. of the eq. (2) formally match for K=3 and the parameters e and 
0o finally read: e = 2 ai/(ao — ai) and O = — 7r )/2- In passing we also quote the major 
semiaxis of the ellipse a = l/V^o — a\. Results of HLSVD are summarized in Tab. 1; for the 
sample test the average estimated values for e and O are respectively (0.93 ± 0.10) x 10~ 3 
and 82° ±5°. 

C. Tagging pixels and radial integration 

We shall focus on the pixel tagging problem. We do not apply data-reduction procedures 
(detector response, diffraction corrections, Bragg peaks removal, etc.) at this stage since 
we are not interested in them. Nonetheless to say that their inclusion will not affect our 
conclusions (see Hammersley et al. 1996 for instance). 

Tagging pixels onto the 2-D detector is a crucial step in the XRPD intensity profiling prob- 
lem. Infact the final XRPD intensity profile follows the radial integration over the ellipses 
domains. Pixels are actually labeled by an adaptive binning algorithm. Adaptive binning is 
an analysis technique that involves subdividing a signal into samples that are more homo- 
geneous than the whole signal (see Sanders et al. 2001 and references therein for instance). 
This technique reveals information about the structure of the signal. Here we have adopted 
the following scheme. The XRPD 2-D image is divided into a set of annuli with common 
centre and eccentricity/phase computed in the previous steps. Then each annulus is tested 
to see whether it meets some criterion of homogeneity {e.g., whether all the pixels in the 
annulus are within a specific dynamic range); here we have chosen a conservative criterion: 
the same statistical representation within each annulus, i.e. the same number of pixels per 
annulus. More complicated criteria can be devised to account for other features (detector 
response and diffraction correction for instance). If an annulus meets the criterion, it is 
not divided any further. If it does not meet the criterion, it is subdivided again and the 
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test criterion is applied to the new annuli. This process is repeated iteratively until each 
annulus meets the criterion. The result might have annuli of several different sizes. Finally 
the 2d scanning is achieved by integrating the 2-D XRPD intensity profile over the annular 
domains tagged in the previous step. The resulting intensity profile corresponding to the 
bottom right image of Fig. 1 is shown in Fig. 2. 

III. RESULTS AND CONCLUSIONS 

We have considered a new approach to the 2$ scanning in powder diffraction experiments 
using 2-D detectors. Our method relies on the decoupling of two main problems: the beam 
centre positioning and the calculation of eccentricity and rotation of the set of ellipses 
recorded. The two tasks are respectively accomplished with by a Simulated Annealing and 
a Hankel Lanczos Singular Value Decomposition, two well known and feasible techniques 
widely applied in several optimization problems. Tagging pixels task is achieved by an 
adaptive binning algorithm by which the whole image is partitioned in a set of annuli with 
the same statistics. Integration over the annular domains finally provides the 2-d intensity 
profile. This approach has been succesfully applied on a set of XRPD CCD images of samples 
of the standard NIST Si 640c collected with a KappaCCD Nonius diffractometer. 
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FIGURES 

FIG. 1. Top left. Original XRPD pattern: the beam centres before and after SA are spotted 
on the image. Top right. Radial distribution function before SA. Bottom left. Radial distribution 
function after SA. Bottom right. Original XRPD pattern: the beam centre after SA together with 
few ellipses tagged by the HLSVD procedure are enlightened. 
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TABLES 

TABLE I. Parameters estimation by HLSVD (K=3) for few ellipses. 



parameters 


ao x 10 4 


a x x 10 s 


fi 


e x 10 2 


<Po 


a 


ellipse 1 


0.117 


0.268 


319° 


0.091 


69.6° 


293 


ellipse 2 


0.119 


0.247 


342° 


0.083 


80.1° 


290 


ellipse 3 


0.120 


0.291 


346° 


0.097 


82.8° 


289 


ellipse 4 


0.146 


0.395 


357° 


0.108 


88.6° 


262 


ellipse 5 


0.145 


0.325 


346° 


0.089 


82.9° 


263 



FIG. 2. XRPD intensity profile corresponding to the bottom right image of Fig. 1. 
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